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CHAPTER  ONE 

INTRODUCTION 

The  need  to  solve  boundary  value  problems  of 
mathematical  physics  has  always  been  appreciated.  From 
the  outset,  enough  assumptions  were  made  to  arrive  at 
a solution  which  approximated  the  physical  behavior  of 
the  problem  in  question.  As  the  geometry  of  the  problem 
became  complex,  approximations  in  the  solution  procedure 
were  accepted  to  arrive  at  a final  result,  leading  to 
graphical  and  numerical  solution  techniques.  Irreconcil- 
ably, only  certain  classes  of  problems  were  admitted 
within  these  techniques,  which  were  basically  determined 
by  the  assumptions  alluded  to  at  the  outset.  The  goal 
of  approximating  physical  behavior  of  this  class  of 
problems  more  closely  or  of  any  problem  not  within  the 
class  can  only  be  realized  by  a reexamination  of  the 
underlying  assumptions,  leading  to  more  complex  governing 
differential  equations.  It  follows  that  the  solution  of 
these  governing  equations  when  considered  over  any  geometry 
other  than  the  most  rudimentary  becomes  cumbersome  and 
even  unobtainable  analytically,  while  the  graphical 
approach  becomes  eliminated  altogether,  leaving  the 
numerical  technique  as  the  apparent  means  of  solution. 

The  numerical  technique  is  flexible  since  it  is 
applicable  to  systems  of  varying  properties  and  nonuniform 
boundary  conditions  besides  being  adaptable  to  complex 
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geometry.  While  the  numerical  technique  involves 
complicated  and  tedious  manipulations,  high-speed, 
large-capacity  digital  computers  have  been  developed  to 
accomodate  them.  Although  numerical  methods  do  not 
provide  general  solutions,  this  could  be  of  little 
consequence,  because  it  may  be  the  particular  solution 
which  is  of  practical  interest.  In  fact,  when  the  general 
solution  is  available,  it  might  prove  to  be  difficult 
and  tedious  to  translate  it  to  a particular  solution. 

Numerical  techniques  have  a rich  history  and  have 
been  extensively  and  successfully  applied  to  the  title 
boundary  value  problems  in  the  form  of  finite  difference 
methods  (FDM)  and  finite  element  methods  (FEM) . There 
are  also  practical  difficulties  associated  with  these 
methods  which  are  also  well  known  [1] . The  fact  that 
we  are  dealing  with  elliptic  partial  differential 
equations  deserves  special  consideration.  It  is  this 
fact  in  combination  with  natural  constraints  on 
physical  constants  that  allows  us  to  relate  interior 
variables  to  boundary  variables,  which  is  the  essence 
of  the  BIE  method.  This  means  that  only  the  boundary  of 
the  problem  need  be  dealt  with  --  a reduction  of  the 
dimensions  of  the  problem  by  one.  Furthermore,  only 
interior  variables  of  interest  need  be  evaluated  if 
desired.  In  order  to  effect  the  BIE  method,  a 
fundamental  singular  solution  of  the  governing  differ- 
ential equations  must  be  found,  which  when  combined  with 


the  desired  solution  through  a reciprocal  relation  that 
will  be  derived  leads  to  integral  equations  defined  on 
the  boundary  of  the  problem  domain.  While  this  represents 
an  exact  means  of  solving  the  system,  it  is  adapted  to  an 
approximate  numerical  procedure  by  a sufficient  number  of 
approximations  within  the  integral  equations  to  convert 
them  to  ordinary  equations  in  unknown  boundary  data. 

These  equations  contain  integrals  of  relevant  variables 
in  which  the  integrals  are  generally  attacked  numerically. 

This  thesis  will  deal  with  the  application  of  the 
BIE  method  to  systems  of  second  order  elliptic  linear 
partial  differential  equations. 

In  Chapter  Two,  we  formulate  the  BIE  for  the  title 
problems  and  indicate  both  the  general  numerical  procedure 
and  specifically  a polynomial  shape  function  (PSFBIE) 
numerical  procedure.  The  details  of  relevant  limit 
processes  needed  for  the  establishment  of  the  BIE  are  given 
in  Appendix  I.  The  details  of  the  calculation  of  the 
matrix  of  coefficients  used  in  the  solution  for  unknown 
boundary  data  for  a polynomial  shape  function  BIE  (PSFBIE) 
with  Legendre  Gauss  quadrature  of  relevant  integrals  is 
given  in  Appendix  II.  Whereas  these  items  are  appendi- 
cized  to  allow  formula tive  ideas  to  flow,  they  constitute 
the  basis  for  actual  numerical  solutions. 

In  Chapter  Three,  the  method  is  applied  to  two 
specific  second  order  linear  elliptic  systems  - anisotropic 
heat  conduction  and  elasticity.  Problems  of  physical 
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significance  in  two -dimens ions  are  solved  via  the  BIE 
method  and  are  compared  with  analytical  solutions,  when 
available,  to  illustrate  the  features  and  reliability  of 
the  method. 

In  Chapter  Four,  we  offer  some  discussion  and 
conclusions  regarding  the  method  in  general  and  some  of 
the  approximations  made  to  convert  the  BIE  to  a numerical 
procedure . 


CHAPTER  TWO 


BOUNDARY  INTEGRAL  EQUATION  FOR 
BOUNDARY  VALUE  PROBLEMS  GOVERNED 
BY  SECOND  ORDER  ELLIPTIC  SYSTEMS 

2.1  ELLIPTIC  PARTIAL  DIFFERENTIAL  EQUATIONS 

2.1.1  Boundary  Value  Problems 

The  theoretical  development  of  the  method  has  been 
given  in  general  form  by  Clements  and  Rizzo  [2] . The 
first  three  sections  of  this  chapter  are  essentially  a 
reproduction  of  material  presented  in  their  paper. 

Consider  the  problem  of  determining  a system  of 
functions  $k(Xj),  (k=l,2,...N,  j=l,2)  throughout  a two- 
dimensional  region  R where  <f>k  satisfies 


3 2 4> 


a. 


ijkil  3Xj3x£ 


— = 0 


i,  k = 1,  2 , . . .N 

j,  l - 1,  2 


Xj  € R 


(2.1) 


and  R is  a simply  or  multiply  connected  domain  with  boundary 
C.  (The  summation  convention  is  implied  on  all  lower  case 
Latin  indices.) 

We  impose  further  that  the  a^^  in  (2.1)  are  real 
constants  satisfying  the  symmetry  condition 


lijkS, 


*klij 


(2.2) 
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The  ellipticity  requirement  can  be  expressed  by 


3*i  3*k 

aijk*  5^7  ^ > ° 


(2.3) 


for  arbitrary  non-zero  3<J)^/3Xj,  3<|>k/3x£. 

We  wish  to  find  a solution  to  (2.1)  which  is  valid 
in  the  two  dimensional  region  R with  boundary  C.  The 
problem  is  posed  with  the  dependent  variables  <}>k  or 
specified  on  C,  where 


3<*>k 

pi  = aijk£  nj 


(2.4) 


in  which  n^  is  the  unit  outward  normal  to  C.  In  other 
words  the  boundary  conditions  can  be  either  of  the 
form  (Dirichlet  type) 


*k(Xj)  = fk(Xj) 


or  of  the  form  (Neuman  type) 


(Xj)  6 C 


Pi(xj)  = gi(x.) 


(x^  e c 


or  of  the  form  (Mixed  problem) 
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^(Xj)  = fk(Xj)  (Xj)  € C1 

Pi(Xj)  = gi(Xj)  (x j ) € c2 

where  C = U c2 • The  functions  f^fxj)  and  g^(xj)  are 
prescribed  on  the  boundary  C,  and  in  accordance  with  a well 
posed  boundary  value  problem,  fv (x . ) and  g. (x. ) are  not 
simultaneously  prescribed  over  the  same  part  of  the  boundary. 

If  the  problem  is  of  the  second  type,  (P^(Xj) 
specified  on  C) , we  require 

' 

P.jdS  * 0 (2.5) 

C 

which  represents  a steady-state  restriction  in  the  scalar 
case  of  N=1  and  static  equilibrium  for  N>1. 

2.1.2  Solution  In  Terms  of  Analytic  Functions 

An  analytic  solution,  by  its  nature,  allows  us  to 
take  derivatives  of  the  fundamental  variable  <t>k  to 
determine  Pi  via  the  constitutive  relation  (2.3)  - making 
the  boundary  integral  equation  formulation  more  concise. 
Accordingly,  we  may  take  the  form  of  the  solution  of  (2.1) 
as 


♦k  = Akf<z> 


(2.6) 
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in  which  f(z),  z = x^  + px2»  is  an  analytic  function  of 
the  variable  z to  be  determined  while  A^  and  p are 
constants.  Substitution  of  (2.6)  into  (2.1)  yields  the 
homogenous  algebraic  system  of  N equations 

(ailkl  + ailk2p  + ai2kl  + ai2k2p2)Ak  = 0 (2'7) 

which  presents  a characteristic  value  problem  if  the 
solution  is  to  be  non-trivial.  Accordingly,  this  requires 
that 

ailkl  + ailk2p  + ai2klp  + ai2k2p2 

The  determinant  in  (2.8)  gives  rise  to  a polynomial  of 

degree  of  2N  in  p,  which  has  N roots,  which  will  be  required 

to  be  distinct.  These  N roots  are  necessarily  complex, 

in  order  that  they  satisfy  the  ellipticity  condition  (2.3). 

Eshelby  et  al  [3]  prove  that  the  ellipticity  condition  is 

violated  for  a real  root  for  the  case  of  N»3.  Furthermore, 

because  the  a.  „ are  real  it  follows  that  the  roots  will 
ijki 

occur  in  complex  conjugate  pairs. 

Denoting  those  roots  with  positive  imaginary  parts 
by  Pa>  a = 1»2,...N,  and  their  conjugates  by  pa,  the 
corresponding  A will  be  identified  as  Aj^  and  A^, 
respectively.  The  real  solution  of  (2.1)  can  then  be 


= 0 (2.8) 


written  as 
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2 Akafa*zo^  + ^ Afcafa(za) 

a=l  a=l 


(2.9) 


where  za  - xi  + Pax2* 

2.1.3  A Reciprocal  Relation  Between  Boundary  and  Interior 
Values 

The  reciprocal  relation  which  we  develop  here  results 
from  an  application  of  Green's  Theorem  to  the  governing 
differential  equation  (2.1)  which  we  will  later  use  in 
conjunction  with  a fundamental  singular  solution  4>k  to 
develop  an  interior  identity.  The  application  of  Green's 
Theorem  is  accomplished  by  letting  4>k  be  a solution  of 


aijk£  3x  3xa  " hi(xl'  x2> 


(2.10) 


with  corresponding  Pk  as  defined  by  (2.4).  Also  let 
4>k  be  a solution  to  a similar  system  with  h^  replaced  by 
h£.  By  the  definition  of  Pk  and  the  divergence  theorem, 
we  may  write 


P^fds 

4 

c 


aijkt  nj*lds 

t 

c 


a-  n 

i;jk£ 


(2.11) 


R 
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Using  (2.10),  this  becomes 


* 

3*  3<tr 

p^rds  = 

ai jki.  3x^  ^dR  + 

h-  d>  r. 


dR 


(2.12) 


By  the  same  operations, 


Pi^dS  = 


3<^k  3<^i 

aijk*  3^  + 


h^dR 


(2.13) 


We  may  now  combine  (2.12)  and  (2.13)  to  obtain 
the  reciprocal  relation 


[Pi*i  “ P^ildS  = 


J 

C 


[h ±<t>r  - h^JdR  (2.14) 


since 


3<J»k  3<^i  3<^i 

ai  j kit  3x£  &x  j aijk£.  3x^  3xj  (2.15 

by  the  required  symmetry  condition  (2.2). 

2.1.4  A Fundamental  Singular  Solution 

The  heart  of  the  BIE  lies  in  finding  a solution 
hi  in  (2.10)  which  satisfies  (2.1)  at  all  points  except 
x * xQ,  (x0€  R)  , whereby  the  singularity  of  the  proposed 
solution  yields  the  necessary  information  about  <P^(xQ) 


to  make  (2.14)  an  interior  identity.  In  other  words,  we 
seek  the  solution  to 


aijk£  3Xj3x£  = Ki6(2  " ?o} 


(2.16) 


x = (xx,  x2) 


xQ  = (a,b) 


where  6 is  the  Dirac  delta  function  and  are  constants 
to  be  determined.  A candidate  for  a solution  would  be 
the  function  of  (2.9)  with 


f (z  ) — ■s — i-  D loq  (z  — c ) 
a v cr  2iri  a A y v a or 


(2.17) 


where  ca  = a + p^b  and  are  constants.  Accordingly, 
(2.9)  becomes 


*k  = 2il  {2AkaDalog(za-ca)  + 2AkaDalog(za-ca)> 


(2.18) 


The  multivalue  properties  of  the  logarithm  provide  that 
transversing  any  closed  path  encircling  the  point  xQ 
causes  ^ to  jump  by  an  amount 


kk  * 2(AkotDa  + AkaDa^ 


(2.19) 
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To  make  identically  zero,  let 


D — 4i  N ■ d . 
a 2 aj  j 


where  d^  are  real  constants  and 


AkaNaj  5kj 


(2.20) 


(2.21) 


The  existence  of  NQj  is  guaranteed  by  the  fact  that 
Aka  is  non-singular.  The  non-singular  character  of  Aka 
is  established  for  the  case  when  (2.8)  has  N distinct 
roots  by  an  extension  of  a result  obtained  by  btroh  [4] 
with  N=3.  By  this  choice  of  Da,  bk  is  identically  zero, 
making  (2.18)  a possible  solution  to  (2.16). 

This  solution  is  valid  everywhere  with  logarithmic 
singularity  at  xQ,  which  is  admissible  within  the  Dirac 
delta  function.  In  order  that  (2.18)  satisfy  (2.16), 
we  must  relate  the  to  the  Da  by  the  residue  of  the 
•integral  over  an  arbitrary  region  enclosing  the  point  xQ. 

We  accomplish  this  by  determining  directly  from 
integration  of  (2.16)  with  ^ as  given  by  (2.18).  By 
this  substitution  (2.16)  becomes 
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2uT  I ailkl  ' ? 


A D 
ka  ci 


a (za-ca) 


A,  D 

_+  2 -JgJl  . 

2 a (Za“ca)2 


ilk2 


_ ^ ^ko^^ci 
_ a (za-ca)J 


i2kl 


- 2 


a.  D 

ka  a 


a (z0-ca) 


+ y AkaDapa 

q (z„-c„)2 


A.  D 

+ ^ ka  a 
a (*0-c0) 


i2k2 


A D p 
y ka  a a 

a ( z- c ) 
L a a 


u - AkaDaPa 
2 « (2a-ca)2 


(2.22) 


Define 


L.  . = (a..,,  + p a . . , _ ) A, 

ija  ijkl  a ijk2  ka 


(2.23) 


From  (2.7) 


(ailkl  + ailk2pa)Aka  (ai2klpa  + ai2k2pa)Aka 


Since  Aka  is  non-zero,  it  follows  that 


L.,  = -p  L . 

lla  a i2a 


(2.24) 


We  can  then  write  (2.22)  as 


L.  D p 

y + y i2ot  ° a \ (2.25) 

111  a (z„-c„)2  a (z„-c„)2  J 
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Consideration  of  equation  (2.16)  shows  that  the  only 
non-trivial  information  determined  by  the  integration  of 
(2.24)  will  be  given  by  the  point  xQ.  Therefore  we  can 
convert  the  region  integral  to  an  arbitrary  contour 
integral  surrounding  xQ.  For  convenience  we  choose  a 
square  of  sides  of  length  two  (with  centroid  ca  which  we 
take  as  the  origin)  as  our  contour.  Green's  theorem 
yields  the  following 


lijk l SXjSx^  dxldx2 


L.  D P 
i2ot  a ct 

zct“ca 


„ Li2aDaP«  1 . 

I *a-Ca'  ^ 


Li2aDa 


z -c 
a a 


Li2aDa 

2.  -= — =~  dx. 
a z -c  -L 

a a 


(2.26) 


The  result  of  (2.26)  will  be  the  following  expression 
for  our  arbitrary  constant  : 


2 L.0  D +2  L. _ D 

^ i2a  a ^ i2a  a 


(2.27) 


or  equivalently 


Ki  * “3d-  2 (Li2aNaj  ” Li2aNaj^dj 
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Thus  (2.18)  is  the  sought  solution  if  the  DQ 
satisfy  (2.20)  and  the  (2.27). 

2.2  THE  BOUNDARY  INTEGRAL  FORMULATION 
2.2.1  General  Formulation  For  Second  Order  Elliptic 
Systems 

We  may  now  use  the  fundamental  singular  solution 
and  the  reciprocal  relation  of  section  2.1.3  to  formulate 
the  boundary  integral  equation.  First  we  note  that  the 
solution  which  we  found  (2.16)  is  still  arbitrary  in 
K..  If  we  let  K.  = <5.  . F,  F an  arbitrary  constant,  in 

i l i;j  2 

(2.27)  for  j=l,2,...N,  we  introduce  a dependence  of 
on  j . Denote  this  dependent  variable  , as  it  is 
obtained  from  (2.18),  and  denote  the  corresponding  Pk 
obtained  from  (2.4)  in  its  dependent  form  as  r^j . Consider 
the  reciprocal  relation  (2.14).  If  we  associate  <f>^  with 
*ij,  P^  with  h£  with  K^6 (x-xQ)  and  let  the  unprimed 

variables  be  a regular  solution  to  (2.1)  (i.e. , h^=0) , 

we  obtain 


[P . ♦ . . - r.  .<fr.  ]ds  - -F 

11]  1]  1 


(x-xQ)dR 


(2.28) 


We  can  realize  an  interior  identity  from  (2.28) 
directly  by  familiar  interpretations  of  the  Dirac  delta 
function.  Alternatively  we  may  realize  it  by  the  classical 


neighborhood  of  the  singular  point  (Figure  1) , and  taking 
the  limits  (see  Appendix  1.1  for  details).  By  either 


FIGURE  I 


4 


17 

We  now  proceed  to  develop  the  boundary  integral 
equation  by  letting  xQ  become  an  element  of  C,  which 
we  will  recognize  xQ  by  our  notation  convention  as  XQ. 
Augmenting  the  contour  by  a rectangle  (Figure  2)  to 
facilitate  easier  integration  and  taking  the  appropriate 
limits  (see  Appendix  1.2)  we  arrive  at  the  boundary  integral 
equation 

r 

[Pi<X)*  (X,  Xo)  - r (X,  Xo)<fr.(X)]ds(X) 

c 

= -X. . (X  ) F$ . (X  ) (2.30) 

i]  ~o  i ~o 


where  X^j  is  related  to  the  inner  angle  0(XQ).  By 
considering  a uniform  distribution  of  <J>^,  which  for 
convenience  we  take  as  unity,  we  may  express  ^ij(xQ)F  as 


X. . (X  ) F = 
1]  -o 


r.  . (X,  X ) ds (X) 
1]  - ~o 


C 


(2.31) 


Substituting  (2.31)  in  (2.30)  we  have  the  boundary 
integral  equation  in  its  most  general  form 


^i  (~o)  ri j (~ ' Vds(?> 


[P.  (X)*..(X,  X ) - r..(X,  X ) <J> . (X)  ] ds  (X) 
1 ~ 13  ~ ~o  13  - ~o  1 - 


J 

c 
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2.3  A NUMERICAL  PROCEDURE 

2.3.1  Numerical  Implementation  of  The  Boundary  Integral 

Equation 

We  note  forthright  that  (2.30)  constitutes  a 
singular  integral  equation  from  which  an  exact  analytical 
solution  could  be  determined  in  principle.  However, 
this  approach  wouxd  be  successful  only  for  the  most 
elementary  problems.  Moreover,  in  such  cases  the  process 
would  still  require  a tremendous  effort  as  compared  to 
other  means  of  arriving  at  the  same  analytical  solution. 
Therefore  the  BIE  is  applied  by  making  a sufficient  number 
of  approximations  to  obtain  a numerically  realizable 
procedure.  Specifically,  the  approach  is  to  approximate 
the  BIE  by  a system  of  linear  algebraic  equations  that 
will  be  solved  numerically  to  obtain  the  desired  values 
of  <f> . (X)  and/or  P (X)  according  to  the  boundary  value 

l ~ i - 

problem  posed. 

As  a means  to  this  end,  four  major  steps  are  taken, 
all  of  which  bear  significantly  on  the  final  result.  They 
can  be  categorized  in  occuring  order  as: 

i)  Discretization  of  the  boundary; 

ii)  Boundary  geometry  representation; 

iii)  Boundary  variable  representation; 

iv)  Evaluation  of  integrals. 

The  discretization  of  the  boundary  is  dividing 
of  the  boundary  curve  into  m segments,  Ca,  over  which 
the  BIE  method  is  applied  (Figure  3 depicts  a typical 
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discretization  pattern).  The  BIE  (2.30)  may  then  be 
written  as 


m 


-A 


ij<?o>F*i<?o>  = 2, 


[pi(X)*ij(X,  XQ) 


- r.  . (X,  X )4>.  (x) ] ds (x) 

ij  - ~o  i - 


(2.32) 


Although  this  step  introduces  no  approximations,  it  has 
a very  significant  effect  on  the  error  which  will  arise 
from  later  approximations,  as  will  be  demonstrated.  There 
is  no  mechanical  procedure  for  this  process,  but  descriptions 
of  the  other  steps  outline  an  "ad  hoc"  set  of  rules  to 
follow  in  discretizing  the  boundary. 

Usually,  no  approximations  of  boundary  geometry 
need  be  made  for  two  dimensional  problems  of  interest; 
however,  to  establish  a solution  method  which  will 
numerically  solve  the  generic  problem,  it  will  be  found 
quite  convenient  to  approximate  the  boundary  geometry. 
Polynomial  shape  function  approximations  are  the  most 
straightforward  approach,  yet  provide  extremely  good 
results  using  low  order  polynomials.  Note  that  the  error 
introduced  here  will  depend  heavily  on  the  boundary 
discretization . 

The  representation  of  boundary  variables  is  the 


first  "true"  approximation  that  need  be  made.  We  require 
it  so  that  we  may  approximately  replace  the  singular 
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integral  equations  involved  with  regular  integral  equations, 
in  order  that  we  obtain  a simultaneous  linear  algebraic 
system  in  the  boundary  variables.  Again,  polynomial 
shape  functions  of  a low  order  lend  themselves  to  very 
good  results  in  this  approximation.  Furthermore,  consider- 
ations of  the  boundary  variables  should  be  made  in  conjunc- 
tion with  considerations  of  the  boundary  geometry  when 
deciding  upon  a discretization  pattern  as  we  note  that 
<t>^  is  a continuous  function,  whereas  is  not  necessarily 
continuous. 

We  now  have  a system  of  regular  integral  equations. 

At  this  point  we  have  a choice  regarding  the  manner  in 
which  we  will  evaluate  these  remaining  integrals.  We  may 
evaluate  them  either  analytically  (employing  the  polynomial 
approximations  made  so  far)  or  by  numerical  integration  of 
a polynomial  interpolating  function.  One  can  certainly 
appreciate  the  complexity  of  the  integrals  involved  (which 
could  very  possibly  lead  to  unobtainable  analytical 
results)  that  is  further  increased  by  higher  polynomial 
approximations  of  boundary  variables,  all  of  which  make 
numerical  integration  more  appealing.  The  choice  of  the 
numerical  integration  procedure  must  also  be  considered 
in  the  discretization  process. 

The  execution  of  these  steps  allow  us  to  arrive 
at  a system  of  n simultaneous  linear  algebraic  equations 
in  n unknowns,  which  can  be  solved  to  obtain  the  desired 
boundary  value  data.  The  n referred  to  here  is  determined 


by  the  number  of  segments  (m)  and  the  degree  of  the 
polynomial  (y)  used  to  approximate  the  boundary  variables. 
2.3.2  Boundary  Integral  Equation  Using  Polynomial  Shape 

Functions  (PSFBIE) 

Inherently  the  BIE  requires  integration  along  the 
contour  (which  we  will  approximate  in  practical  problems 
of  interest)  that  inevitably  leads  to  difficulties  in 
integration.  Consequently,  in  the  process  of  approximating 
boundary  geometry  and  boundary  /ariables  there  is  a distinct 
advantage  in  changing  the  variable  of  integration  by  mapping 
each  segment  to  a straight  line. 

This  mapping  is  accomplished  by  expressing  the 
cartesian  coordinates  of  each  point  on  the  given  segment 
in  terms  of  the  cartesian  coordinates  of  the  nodes  and 
the  coordinate  £ in  the  mapped  space  [1] . Similarly  we 
express  the  boundary  variables  of  each  point  in  terms  of 
the  boundary  variables  of  the  nodes  and  the  coordinate  £ 
in  the  mapped  space.  Mathematically  we  may  express  this  as 


Xj  (Q) 

= Ma(5)Xj 

Pi  (Q) 

- Ma(5)P? 

j = 

*i(Q) 

= Ma(£)<J»2 

i — 1,2, ...N 

(2.33) 

where  Q is  a point  on  the  boundary  where  we  are  integrating, 
Ma  is  the  polynomial  of  degree  y which  is  used  to  approxi- 
mate the  said  variables,  and  a has  the  range  of  the  number 
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of  nodes  on  the  segment  (y  + 1) . Figure  4 shows  a 
graphical  interpretation  of  the  first  relation  of  (2.33). 
The  coefficients  of  the  polynomials  Ma  depend  upon  the 
locations  of  the  nodes  in  the  mapped  space.  Logically, 
we  would  locate  the  nodes  symmetrically  on  the  straight 
line  in  tne  mapped  space.  Furthermore,  we  will  define 
5=0  as  the  midpoint  of  the  straight  line. 

a 

5—1  5=0  5=1 

X| 

FIGURE  4 

The  mapping  procedure  outlined  so  far  has  an 
associated  Jacobian  which  is  used  in  the  evaluation  of 
integrals  in  the  real  space.  In  other  words,  ds(X)  must 
be  expressed  as  J(5)d5*  To  obtain  J(5)  consider  the 
differential  elements  depicted  in  Figure  5 which  are 
related  as 

ds2  = dx£  + dx| 


(2.34) 


FIGURE  5 


Using  (2.33),  the  differential  dX . is 

3 


a 


_ dM  ( £ ) 


dXj  = — d r1- 


for  which  the  case  of  j=l,2  yields 


ds  = 


= J(5)dC 


(2 


Equation  (2.32)  may  now  be  expressed  as 


-X 


m 


ij<?o>F+i<5o>  * 2 
0=1 


{Ma(C)  [P“a*i;j  (X(£)  , XQ) 


- r._.(xu),  x0)$“a]}j(5)dc  (2 
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where  is  interpreted  as  the  i component  of  P on 

segment  o at  local  node  a and  <)>^CT  has  a similar  inter- 
pretation. 

2.3.3  Numerical  Quadrature  of  Integrals  in  the  PSFBIE 

A broad  definition  of  numerical  quadrature  may  be 
taken  as  the  numerical  integration  of  a polynomial  inter- 
polating function,  which  will  yield  an  exact  result  if  the 
function  to  be  integrated  is  of  the  same  order  as  the 
interpolating  function  or  less.  Mathematically,  this 
can  be  expressed  as 

» 

f(5)dC  = [ 2 w f(?  )]AC  (2.38) 

i=l  1 

AC 

where  AC  is  the  interval  over  which  we  are  integrating, 
v is  the  number  of  function  evaluations  in  AC,  and 
Cj_  are  the  roots  of  the  interpolating  function  (in  the 
interval  in  which  it  is  valid)  with  appropriate  weights 

wi* 

For  quadrature  formulae  of  the  form  (2.38)  the  most 
accurate  are  Gaussian  quadrature  formulas  [5] . The 
particular  type  of  Gaussian  quadrature  is  determined  by 
the  interpolating  polynomial  associated  with  it.  As  yet, 
the  length  of  the  straight  line  in  section  2.3.3  purposely 
has  not  been  specified  with  the  foresight  of  choosing  it 
to  obtain  the  most  accurate  numerical  quadrature  in  the 
mapped  space.  Accordingly,  the  choice  of  Legendre -Gauss 
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Quadrature  determines  the  straight  line  as  the  interval 
£=-1  to  5*1 • This  choice  leads  to  the  exact  integration 
of  polynomials  of  degree  2v-l  or  less  where  v is  the 
number  of  function  evaluations  in  (2.38). 

Note  that  (2.35)  contains  singular  integrals  when 
XQ  is  on  segment  a.  In  this  case  Legendre-Gauss  Quadrature 
will  not  give  accurate  results.  Consequently,  modifications 
must  be  made,  for  which  the  reader  is  referred  to 
Appendix  II.  Appendix  II  also  describes  the  evaluation 
of  non-singular  integrals  involved  in  the  calculation  of 
the  coefficient  matrix. 

2.3.4  Assembly  and  Solution  of  the  Coefficient  Matrix 
We  now  seek  to  represent  the  Boundary  Integral 
Equation  in  its  present  form  (2.37)  as  a matrix  equation 
so  that  we  may  solve  it  as  a simultaneous  linear  algebraic 
system.  Note  that  by  using  (2.31),  equation  (2.30)  can 
be  expressed  as 

' 

[<J>  (X)  - <J>  (X  ) 1 r (X,  x ) ds (x) 

1 ~ 1 ~o  1 j ~ ~o 

C 

' 

* X0)ds(X)  (2.39) 

C 


whereby  its  counterpart  in  the  PSFBIE  equation  (2.37) 


assumes  the  form 


28 


m 

2 < 

0=1 


*“aMa(£)r..(X(£),  Xo)J(C)dC 


Lcr 


- *i<*o> 


rij(x(C) , xo)J(C)d 5 


m 

2 

0=1 


.ao..a 


pi  M (X(C)  , XQ)  J(C)dC  (2.40) 


When  calculating  the  indicated  integrals  (either  analytically 
or  numerically),  equation  (2.40)  can  be  concisely  formulated 
in  matrix  form  as 

2 - »i(X0)PX?j]  = 1 P?°b«[  (2.41) 

0=1  J a=l 


or 


[A]  { <J) } = [B]  {P} 


(2.42) 


where  [A] 

{*} 

[B] 

{P> 


'• “ - «I3i 

[b“i 

{p“0} 
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(2.43) 

(2.44) 

(2.45) 

j| 

The  details  of  the  calculation  of  a£ j , j and  b£ j for 
both  singular  (i.e.,  XQ  on  segment  a)  and  non-singular 
integrals  evaluated  numerically  are  supplied  in  Appendix  II. 

The  elliptic  system  (2.1)  will  then  be  solved  according 
to  the  form  of  the  boundary  conditions. 

For  the  Dirchlet  problem,  (2.42)  becomes 

{c}  = t B ] {P}  (2.46) 

where  {c}  = [A]  {<{>}. 

For  the  Neuman  problem  we  have 

[AH*}  “ {c}  (2.47) 

where  {c}  = [B]{P>.  However,  Neuman  problems  do  not  have 
unique  solutions,  making  the  [A]  matrix  singular.  To 
deal  with  this  situation,  we  constrain  one  value  of 

arbitrarily,  so  that  we  may  uniquely  obtain  all 
other  values  of  {<j>?°}.  This  could  be  realized  by 
specifying  one  {<fr®a}  to  be  zero  and  then  eliminating 


’ 

a“°  = Ma(5)r  (X(C),  XQ)J(S)d5 

J Jc  J 

a 

FXij  = rij (X(5) » XQ)J(C)dC 

4 C 

a 

b*a.  = Ma(£)$  (X(C),  X )J(C)dC 

1]  1]  ~ ~o 

C 

a 
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the  corresponding  row  and  column  in  [A]  and  the  corresponding 
element  of  {c}.  Thus  we  deal  with  the  reduced  system. 

tA*]{<fr*}  = {c*>  (2.48) 


Altneratively , the  value  of  one  element  of  t<j)“a} 
may  be  specified  arbitrarily  and,  considering  the  corres- 
ponding {P?a}  to  be  unknown,  the  system  can  be  solved  as 
a mixed  problem  as  will  be  outlined  next. 

For  a mixed  problem,  (2.42)  can  be  expressed  as 

[A*]  = {c*}  (2.49) 


where  [A*J  = (a? j ] 

{$*}  = U“0*} 

{ c*  } = [B* ] {P* } = [b“°*](p“a*> 


l<i , j<N 


The  starred  quantities  are  defined  per  the  known  boundary 
variables  on  segment  a.  For  a segment  where  is 
defined 
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The  previously  unknown  boundary  data  in  equation 
(2.46),  (2.48)  or  (2.49)  may  now  be  solved  for  by  any  of 
various  standard  methods. 

Once  the  unknown  boundary  data  is  obtained  by 
one  of  these  means,  the  <J>j  (XQ)  in  the  interior  identity 
(2.29)  may  be  evaluated  directly  from  a polynomial  shape 
function  adaptation  of  (2.29). 


CHAPTER  THREE 


TEST  PROBLEMS 

3.1  PLANE  ANISOTROPIC  HEAT  CONDUCTION 

3.1.1  Application  of  Boundary  Integral  Equation  Formulation 
For  the  problem  of  plane  anisotropic  heat  conduction 

the  governing  differential  equation  (2.1)  takes  the  form 


3 


ji. 


= 0 


j,*  = 1,2  (3.1) 


where  our  symmetry  requirement  (2.2)  provides  that 
a_.^=a^_..  The  characteristic  parameter  P^  may  be  solved 
from  (2.8)  as 


-a 


Pa  = 


12 


± r 


l12 


2-  a 


ll“22 


a 


22 


(3.2) 


The  determination  of  A from  (2.7)  is  now  arbitrary, 
which  we  choose  as  unity  for  convenience.  It  then 
follows  that  N is  also  unity  by  (2.21).  The  L of 
(2.23)  are  then 


Lj  " (ajl  * Vj2>  (3'3) 

The  only  Lj  of  significance  is  L2 , which  is  used  for 
the  determination  of  d in  terms  of  K,  viz. 
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d = 


-2K 


a22i^Pa  " Pa> 


(3.4) 


Finally  D is  determined  from  (2.20)  as 


a22(pct  ~ 


(3.5) 


and  the  fundamental  solution  is  assembled  per  equation 
(2.18).  Now  K is  taken  as  unity  per  section  2.2.1  and 
<j>  becomes  associated  with  $,  while 

I 


r = 


3 

lj  l 3x(i 


(3.6) 


to  arrive  at  the  boundary  integral  equation. 

3.1.2  Solved  Problems 

The  purpose  of  this  section  is  twofold.  First  and 
foremost,  it  is  to  demonstrate  the  ability  of  the  method 
to  solve  problems  of  physical  significance  reliably  and 
secondly  to  show  the  limited  amount  of  input  required  in 
both  data  preparation  and  actual  computer  input.  To 
accomplish  this,  two  problems  are  solved:  one  for  which 
the  exact  solution  is  known,  where  the  comparison  can  be 
made  directly;  and  the  other  one  which  by  its  nature 
allows  us  to  draw  on  our  intuition  to  assess  the 
reliability  of  the  solution.  Furthermore,  essential 
parameters  of  the  particular  problems  are  varied  to 
demonstrate  the  ability  of  the  method  to  discern  the 


i 


' 


jk 
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variation. 

(a)  Consider  a linear  temperature  distribution, 
4>=cx1+d,  across  a rectangular  region  of  dimensions 
21  x 2w  as  shown  in  Figure  6.  Now  P can  be  expressed  as 


P 


allcnl  + a12cn2 


Since  we  have  analytical  expressions  for  both  temperature 
and  heat  flux  as  functions  of  position,  we  will  consider 
the  problem  as  a mixed  one. 

The  nature  of  this  problem  also  lends  itself  to 
be  used  to  demonstrate  another  important  aspect  of  the 
solution  procedure  which  heretofore  has  received  little 
attention,  namely  the  accuracy  of  the  method  as  the 
system  becomes  weakly  elliptic.  The  weakening  of  the 
ellipticity  is  apparent  from  equation  (3.2)  as  the 
discriminant  approaches  zero.  Thus  the  approach  is 

strictly  a function  of  the  material  parameters  a . . . 

D * 

In  order  that  we  may  follow  this  approach  analyti- 
cally, we  set  a11  to  unity,  a12  to  zero,  leave  a^2  as 
a variable  and  rotate  the  material  through  an  arbitrary 
angle  9.  The  transformed  material  parameters  are  given 
generally  as 
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an 

aos20 

sin20  2sin0oos0 

all 

a22 

sin20 

cos2  9 -2sin0cos0 

< 

a22 

a12 

V J 

-sin0cos0 

sin0oos0  (oos20-sin20) 

a12 

^ / 

Investigating  da^/^®  to  find  0 such  ti  a12  is  a 
maximum,  which  would  cause  the  system  to  become  parabolic, 
we  find 


da£2 

d0 


0 


sin20 


cos20  + a22 (cos20-sin20) 


which  leads  to 


cos20  = 1/2 

indicating  that  the  maximum  occurs  at  0=ir/4 . Substituting 
this  value  into  (3.7)  leads  to 


n 

rH 
' .-1 

(a22 

+ l)/2 

II 

CM 
' CM 

(a22 

+ l)/2 

12  = 

(a22 

- l)/2 

so  that  in  the  limit  as  a approaches  infinity  all  of 

* 

the  transformed  parameters  become  equal,  thus  indicating 
a parabolic  system. 


We  arbitrarily  set  dimensions  l and  w to  unity  and 
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take  the  constants  of  the  exact  solution  <p,  c and  d,  also 
to  be  unity  while  a22  is  varied  for  0=ir/4 . The  small 
amount  of  input  required  can  be  seen  by  using  the  crudest 
possible  discretization  for  a quadratic  shape  function  BIE 
(QSFBIE)  of  one  segment  per  side  for  a total  of  eight  nodes. 

The  problem  is  posed  with  temperature  specified 
on  x^  = +1  and  flux  specified  on  x2  = +w,  for  which 
boundary  solution  comparisons  are  summarized  in  Table  1 
and  depicted  graphically  in  Figure  7.  All  solutions 
were  obtained  using  a 4-point  Legendre  Gauss  rule  for 
numerical  quadrature.  The  average  absolute  error  (6^) 
of  Table  1 is  defined  as 


BIE 


- ¥ 


exact ‘ 


exact 


(3.8) 


where  V is  either  variable  being  solved  for  in  the  mixed 
problem  (T  or  3T/3n)  and  is  the  number  of  nodes  at 
which  the  particular  variable  is  solved  for. 

(b)  Consider  the  problem  of  a rectangular  region 
subject  to  uniform  temperatures  at  two  opposing  sides  and 
insulated  on  the  other  two  sides  as  depicted  in  Figure  8. 
The  problem  is  a mixed  type  for  which  there  is  no  known 
analytical  solution.  However,  in  the  general  case  of  plane 
anisotropic  heat  conduction,  the  specifics  of  the  problem 
allow  a limited  analysis  which  provides  a "check"  of  sorts 
on  the  solution. 


38 


TABLE  1 

AVERAGE  ABSOLUTE  ERRORS  x 104 
9 = it/ 4 


a 

22 

all 

DISCRIMINANT 

FLUXES 

TEMPERATURES 

2 

2. 

12 

2 

3 

3. 

16 

2 

4 

4. 

18 

7 

5 

5. 

24 

8 

6 

6. 

28 

3 

7 

7. 

33 

7 

8 

8. 

52 

24 

9 

9. 

73 

47 

10 

10. 

93 

76 

11 

11. 

114 

109 

12 

12. 

113 

147 

13 

13. 

174 

188 

14 

14. 

176 

234 

15 


15 


196 


282 


T r 


C/5 


FIGURE  7 
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Consider  the  intersection  of  an  isotherm  with  an  arbitrary 
contour  as  shown  in  Figure  9.  If  we  define  the  intrinsic 
coordinate  s along  the  isotherm,  the  chain  rule  for 
differentiation  can  be  used  to  express  the  direction  of 
the  variation  of  T with  s as 


3T  = 3T  ^Xi  + 3x2 

<Ts  3x^  3s  3x  3s 


(3.9) 


where  the  position  vector  x(s)  of  an  arbitrary  point  (Q) 
with  respect  to  the  origin  is  given  by 


x(s)  = x (sje^  + x2(s)e2 


(3.10) 


and  ej  is  a unit  vector  basis  defined  in  the  direction 
of  the  coordinates  x^ . Along  an  isotherm  3T/3s  is 


FIGURE  9 
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1 


necessarily  zero  and  3x^/3 s will  be  a unit  vector  tangent 
to  the  isotherm  which  we  will  denote  n*.  Thus  equation  (3.9) 
yields 


3T 

3x 


1 


nj 


-3T 

3x2 


nS 


j * 1,2  (3.11) 


The  definition  of  an  isotherm  provides  that  it  is  a line 
of  constant  flux,  the  value  of  which  can  be  expressed  on 
the  boundary  as 


3* 


3T  _ 


(3.12) 


where  k is  a constant  and  n^  are  components  of  the  outward 
normal  of  the  boundary.  Now  we  satisfy  equation  (3.11) 
in  order  to  use  the  relation  (3.12).  Choosing 


(3.13) 


equation  (3.11)  is  satisfied  identically  and  cam  be  used 
to  determine  the  direction  in  which  the  isotherms  intersect 
the  boundary  C. 

We  apply  this  analysis  to  the  problem  in  question 
by  noting  that  the  two  opposing  sides  maintained  at  uniform 


temperatures  will  have  n*  parallel  to  the  outward  normal 
of  the  boundary  while  along  the  insulated  sides,  k in 
equation  (3.12)  is  zero.  Because  k is  zero,  we  may  express 
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equation  (3.12)  along  sides  AB  and  CD  as 


3T  ^ 3T 

a21  3x-i  a22  5xT  “ 0 


Using  equation  (3.13)  we  obtain 


-a21n2  + a22nl 


= 0 


which  we  satisfy  by  choosing 


n*  = a 
1 21 


n*  = a„„ 
2 22 


(3.14) 


(3.15) 


(3.16) 


so  that  the  direction  of  the  directional  derivative  (3.9) 
is  given  by  equation  (3.16).  Note  that  the  direction 
n*  is  the  same  along  AB  and  CD. 

To  demonstrate  the  ability  of  the  method  to 
recognize  variations  in  parameters,  the  problem  is  posed 
similar  to  problem  (a)  (i.e.,  a11=l,  a12=0,  a22=a22' 

A=w=l)  except  that  this  time  a22/all  is  ^xed  and  9 
varied  so  that  equations  (3.16)  can  be  expressed  as 


n*  = a ' - (a22  - a21> cosOsine 

n2  ’ a22  = a22COs!0  + an3in29 


(3.17) 


1 
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The  stability  of  the  solution  procedure  is  demonstrated 
by  choosing  a22  as  10  (representing  a highly  anisotropic 
material) . Note  that  this  problem  as  posed  reduces  to 
problem  (a)  if  0=mir/2.  This  check  can  be  realized  only 
approximately  by  a numerical  procedure;  however,  the  fact 
that  n*  is  the  same  for  a given  isotherm  along  AB  and  CD 
can  be  seen  from  the  graphical  interpretation  of  the 
results  (Figures  10-19) . 

Figures  10-14  are  temperature  contour  maps  of  a 
40  node  boundary  solution  and  subsequent  interior  point 
evaluations  (a  9x9  square  gridwork  of  points  within  the 
geometry  of  the  problem) . Continuity  of  contours  was 
obtained  by  linear  interpolation  of  the  above  data  points 
in  both  the  x^  and  the  x2  directions.  The  values  marked 
within  the  different  regions  designate  the  temperature 
ranges  given  in  Table  2 . 

Figures  15-19  are  isometric  plots  of  temperature 
-•us  position  viewed  from  four  different  azimuths  (i|>) 
foi  clearer  interpretation.  Nodes  of  the  gridwork 
represent  uniformly  spaced  points  at  which  the  temperature 
was  evaluated  numerically  (with  the  exception  of  lines 
AD  and  BC,  where  the  temperature  was  specified) . Nodes 
along  lines  AB  and  CD  are  boundary  solution  points  and 
all  other  nodal  temperatures  were  evaluated  subsequently 
by  the  interior  identity  by  using  the  boundary  solution. 
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TABLE  2 

LEGEND  FOR  FIGURES  10-14 


VALUE 

0 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


TEMPERATURE  RANGE 
0.00-0.15 
0.15-0.35 
0.35-0.55 
0.55-0.75 
0.75-0.95 
0.95-1.15 
1.15-1.35 
1.35-1.55 
1.55-1.75 
1.75-1.95 
1.95-2.00 


FIGURE  12 


D 
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FIGURE  17  d 
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All  connections  between  gridwork  nodes  were  made  by  spline 
function  fits  along  lines  of  constant  xi  and  all  solutions 
were  obtained  from  a QSFBIE  with  a 4 -point  Legendre  Gauss 
rule  for  numerical  quadrature. 


3.2  PLANE  ANISOTROPIC  ELASTICITY 

3.2.1  Application  of  Boundary  Integral  Equation  Formulation 
For  the  problem  of  plane  anisotropic  elasticity 
the  governing  differential  equation  is  of  the  form 


3 2<J>k 

aijk£  3x . 3x  = 

3 S, 


j,*  = 1,2 

i,k  = 1,2, ...N 


where  N can  be  2 or  3.  The  significance  of  the  two 
possible  values  for  N lies  in  the  fact  that  the  principal 
axes  of  the  material  need  not  be  parallel  to  the  x,  -x2 
plane.  In  other  words,  the  three  rotational  angles 
necessary  to  describe  the  material  constants  in  terms  of 
the  problem  coordinates  Xj  can  all  be  varied.  Consequently, 
for  N=3,  the  antiplane  problem  is  also  solved  in  addition 
to  the  plane  problem  (N=2) . We  note  here  that  if  a 
plane  problem  is  posed  as  a problem  with  N=3  there  will 
be  repeated  roots  of  (2.8)  making  the  characteristic 
vectors  AkQ  of  (2.7)  non-unique.  However,  the  plane  and 
antiplane  problems  in  this  case  are  uncoupled  and  two 
of  the  characteristic  vectors  can  be  chosen  arbitrarily. 
Furthermore,  note  that  for  the  case  of  an  isotropic 
material  the  characteristic  roots  pa (o=l , 2 , . . .N)  occur 
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as  repeated  roots  from  (2.8).  Consequently  the  are 
not  unique,  nor  is  there  any  uncoupling,  leading  to  a 
breakdown  of  the  procedure. 

As  in  the  case  of  plane  anisotropic  heat  conduction, 
the  sequence  of  solution  of  the  characteristic  variables 
(Pa/  Ak(J,  Na j , Lija,  dj , and  Da)  is  the  same;  but  the 
complexity  of  the  calculations  do  not  allow  us  to  write 
out  explicit  expressions  for  them.  The  formation  of  the 
kernel  is  then  straightforward  per  (2.18),  but  the 
assembly  of  the  T—  by  (2.4)  involves  four  nested  loops 
which  may  be  reduced  to  two  by  the  following  procedure 
utilizing  L^ja.  We  have  identified  the  real  solution 
4>k  by  (2.9)  as 


*k  ■ 2Val!«)  + 2Aka£a(V 

ct  a 


(3.18) 


Now  can  be  written  from  v 2 . 4 ) as 


(3.19) 


3<J>k  ^k 

Pi  ^aijkl  3*1  + aijk2  3*2^ 


Taking  the  derivatives  indicated  by  (3.19)  on  <J>k  as 
given  by  (3.18)  yields 


Pi  Jt(aijkl  + aijk2pa)Akcx]  fa(za)nj 

+ J[(aijkl  + aijk2pa)^ka1^ci(za)nj  (3,20) 


Hfcirfr  It  i 
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We  recognize  the  quantities  in  brackets  as  L. • so  that 

x j va 

we  may  write  (3.20)  as 


Pi  - 1 2 Lijafa<za)nj  + 2 l"j  (3‘21» 


Finally 


= 


a 


alza'  2ttT  (za  - ca) 


T'T  \ - 1 a 

a Zcr  2iri  (za  - ca) 


so  that  can  be  expressed  as 


if  Da  Da  1 

Pi  = 2¥I  |_  ^ Lija  (Zo-Co)  - 2 Lija  (za-ca)  ] nj 


(3.22) 


Thus  when  we  assemble  I*  — , it  assumes  the  form 


r. 


[2  l, 


a j 


ij  2iri-  |_~  ika  (za-ca) 


^ Lika  3 nk 


(3.23) 


3.2.2  Solved  Problems 


In  this  section  we  wish  to  demonstrate  the  ability 


of  the  method  to  solve  problems  of  physical  significance 
reliably  and  indicate  the  stability  of  the  method  as  the 
governing  differential  equation  we  are  trying  to  solve 
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becomes  parabolic.  As  in  the  problems  solved  in  section 
3.1.2,  essential  parameters  of  the  particular  problems 
are  varied  to  demonstrate  the  ability  of  the  method  to 
discern  the  variation.  The  solved  problems  considered  in 
this  section  deal  only  with  the  case  of  N=2 . 

(a)  Consider  a general  anisotropic  rectangular 
plate  of  dimensions  21  x 2w  in  simple  tension  as  shown 
in  Figure  20.  The  exact  solution  for  Jl=w=a  with 


(-a, -a)  = u2 (-a, -a)  = ui,2 (“a'x2^  = 0 


is 


ui - 


su(xi + 


a) 


= 


( s 12 (x2  + a) 


S16 (xl  + a) * 


where  s_  represent  material  compliances  in  contracted 
notation.  Note  that  for  the  case  of  general  anisotropy 

u2  = u^x^  x2) 

because  of  the  s^g  term.  A non-zero  slg  term  may  be 
obtained  by  rotating  any  material  other  than  isotropic 
(or  transversely  isotropic  about  the  x1~x2  plane)  through 
am  arbitrary  angle  0 about  the  x-j  axis.  In  certain 
material  classes  8=m7r/2  (m=0 , 1 , 2 , . . . ) will  not  provide 
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a non-zero  term,  but  such  cases  will  not  be  considered 
here.  Accordingly,  the  simplest  choice  is  that  of  a 
transversely  isotropic  material  which  does  not  have 
symmetry  about  x2=0.  This  material  will  have  5 independent 
material  constants,  only  four  of  which  will  enter  into 
the  x^-x2  plane  problem,  namely  c^^,  c2222'  c1122'  an<* 
c1212*  For  conven^ence  we  shall  denote  them  in  the 
following  manner 


The  determinant  polynomial  (2.8)  for  this  type  of  material 
aligned  in  the  x1~x2  coordinate  axes  (i.e.,  8=0)  may  be 
written  in  explicit  terms  as 


ALp4  - (F2  + 2FL  - AC) p2  + CL  = 0 (3.24) 


The  biquadratic  roots  of  (3.24)  are 

_2  _ (F2+2FL-AC)  + i/4ACL2- (F2+2FL-AC) 2 , , 

p K3.Hl 


Examination  of  the  approach  of  the  discriminant  of  (3.25) 
to  zero  shows  the  approach  of  the  elliptic  system  of  this 


particular  problem  to  becoming  parabolic.  We  may  observe 
the  approach  analytically  by  putting 

A = n 
c = 3n/2 
f = n/io 
l = n/5 

and  let  , 8 , 7 , . . . 1 to  cause  the  discriminant  to  approach 
zero  while  satisfying  the  ellipticity  requirement 
(equation  (2.3)).  In  such  a manner,  the  accuracy  of 
the  solution  procedure  as  well  as  the  stability  of  it 
can  be  demonstrated. 

The  results  of  this  survey  with  a=a=l  are  summarized 
in  Table  3.  The  average  absolute  error  (£^)  has  a similar 
interpretation  as  in  section  3.1.2  except  that  4*  of 
equation  (3.8)  represents  either  tractions  (t)  or 
displacements  (u)  here.  All  solutions  were  obtained 
using  one  segment  per  side  in  a QSFBIE  (8  nodes  total) 
with  4-point  Legendre  Gauss  rule  for  numerical  quadrature. 

(b)  Consider  an  elastic  anisotropic  plate  bounded 
by  two  hyperbolas  and  two  equal  rectilinear  sections 
subject  to  uniform  tensile  loading  (a)  in  the  direction 
of  the  normal  on  the  rectilinear  sections  (Figure  21) . 

For  a plate  of  infinite  extent  in  the  x^^  directions, 
with  axes  as  shown,  the  hyperbolic  edges  will  be  given 
by 


68 


TABLE  3 

AVERAGE  ABSOLUTE  ERRORS  X 104 


T 

\ 

II 

<r> 

n 

DISCRIMINANT 

TRACTIONS 

DISPLACEMENTS 

9 

-12220. 

122 

116 

8 

-7628.8 

122 

116 

7 

-4471.9 

122 

116 

6 

-2413.8 

122 

116 

5 

-1164.1 

122 

116 

4 

-476.8 

122 

116 

3 

-150.86 

122 

116 

2 

-29.80 

122 

116 

1 

-1.8625 

122 

116 

x2  X2 
X2  _ X1 


The  variable  of  interest  in  this  problem  is  the 
distribution  of  at  the  narrowest  section  (i.e.,  along 
x^=0) . The  geometry  of  the  plate  is  characterized  by 
defining  c=a/b,  which  is  the  governing  parameter  of  the 
stress  distribution  in  the  plate  for  a given  material. 

An  approximate  solution  for  the  stress  distribution 
of  the  general  anisotropic  plate  is  given  by  Lekhnitskii  [6] 
from  which  we  compare  normal  stress  distributions  along 
the  narrowest  cross  section.  If  we  choose  an  orthotropic 
material  with  its  principal  directions  (i.e.,  fiber 
direction  angle  0 of  Figure  21  equal  to  mir/2 , m=0,l,2,...) 
we  may  use  quarter-symmetry  to  pose  the  problem  as 
depicted  in  Figure  22  and  obtain  the  normal  stress 
distribution  directly  from  the  boundary  solution  tractions. 
The  parameter  xj  (Figure  22)  represents  the  finite  length 
of  the  plate  used  for  the  numerical  solution  modelling. 

The  BIE  solution  is  compared  with  the  approximate 
solution  for  an  orthotropic  material  having  stiffnesses 
(in  units  of  lb/  [in2 • 10s  ] ) c11=1.2035,  c22*0 . 60176 , 
c12=0. 04276,  and  c66=0.0699,  which  are  constants  (cf.[6]) 
for  a type  of  plywood.  The  problem  was  posed  with 
a=a=xj=l,  c=l/10,  1/2,  1,  2 for  9*0  and  8=ir/2.  Comparison 
stress  distributions  are  shown  in  Figures  23-26,  in  which 
the  applied  traction  is  shown  on  the  right  hand  side 


FIGURE  21 
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to  serve  as  an  indication  of  the  magnitude  of  the  stress 
concentration.  The  reliability  of  the  method  is  emphasized 
by  the  crude  discretizations  of  the  problem  of  one  and 
two  segments  per  side  in  a QSFBIE  (8  and  16  nodes  respect- 
ively) . All  solutions  were  obtained  using  a 4 -point 
Legendre  Gauss  rule  for  numerical  quadrature. 
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LEKHNITSKII 
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FIGURE  24b 


CHAPTER  FOUR 


DISCUSSION  AND  CONCLUSIONS 

The  purpose  of  the  work  in  this  thesis  is:  (1)  to 
develop  the  algorithm  derived  by  Clements  and  Rizzo 
for  the  solution  of  boundary  value  problems  governed  by 
a system  of  second  order  linear  elliptic  partial  differ- 
ential equations  into  an  explicit  numerical  procedure; 

(2)  verify  the  numerical  procedure  by  comparing  numerical 
solutions  by  it  with  exact  solutions. 

The  first  purpose  is  motivated  by  the  need  to 
solve  the  title  problems  numerically  coupled  with  the  fact 
that  the  boundary  integral  equation  method  is  especially 
well  suited  for  problems  governed  by  linear  elliptic  partial 
differential  equations  [7] , [8] . The  advantage  of  the  BIE 
method  over  other  numerical  methods  such  as  the  Finite 
Difference  Method  (FDM)  and  the  Finite  Element  Method  (FEM) 
in  attacking  the  title  problems  is  that  '.he  numerical 
process  of  the  BIE  method  is  involved  only  with  the  one- 
dimensional bounding  curve  or  curves.  This  serves  to 
greatly  reduce  the  amount  of  data  preparation  and  computer 
core  and  time  requirements,  especially  when  the  ratio  of 
the  circumference  to  the  area  of  the  domain  is  relatively 
low.  The  amount  of  reduction  is  dependent  on  the  degree 
of  the  polynomial  shape  function  (y)  necessary  to 
approximate  boundary  geometry  and  boundary  variables 
appropriately.  In  practice,  quadratic  shape  function 
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r 


approximations  represent  boundary  geometry  and  boundary 
variables  well,  while  providing  adequate  convergence 
for  decreasing  mesh  length  interval  [9] , [10] . The  error 
analysis  for  the  piecewise  quadratic  Lagrange's  polynomial 
approximation  [11]  shows 


I I f - P 


< 


where  f is  the  function  to  be  approximated,  p2  is  the 
piecewise  quadratic  Lagrange's  interpolating  polynomial, 
and  h is  the  mesh  length  interval  between  nodes. 

One  should  note  that  although  computer  core 
requirements  will  be  less  for  the  BIE  method  than  for  a 
numerical  method  such  as  the  FEM,  the  coefficient  matrix 
involved  (cf.  section  2.3.4)  is  a full  matrix  unlike  the 
banded  sparse  matrix  of  the  FEM.  This  is  an  undesirable 
computat iona 1 feature . 

As  to  the  second  purpose,  the  test  problems  of 
Chapter  3 demonstrated  the  ability  of  the  numerical 
procedure  to  solve  problems  of  significance  reliably. 
Moreover,  the  accuracy  of  the  solution  procedure  as  the 
partial  differential  equations  became  weakly  elliptic 
was  demonstrated  through  analytic  surveys.  All  the 
problems  presented  were  run  on  an  IBM  370  Model  165  computer. 

It  is  appropriate  at  this  time  to  make  some  remarks 


regarding  the  algorithm  derived  by  Clements  and  Rizzo. 
The  most  significant  aspect  of  the  algorithm  is  its 
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generality  of  application  to  second  order  linear  elliptic 
systems.  Furthermore,  the  conciseness  of  the  formulation 
makes  the  numerical  interpretation  of  the  algorithm  a 
straightforward  procedure.  However,  the  algorithm  does 
not  represent  the  first  application  of  the  BIE  method  to 
problems  involving  anisotropy. 

The  application  of  the  BIE  to  plane  anisotropic 
elastic  bodies  was  first  performed  by  Rizzo  and  Shippy  [12] . 
Their  formulation  was  based  on  the  stress  field  due  to  a 
point  force  in  an  infinite  sheet  of  anisotropic  elastic 
material  obtained  by  A.E.  Green  [13].  Green's  point  force 
solution  is  derived  from  a stress  function  form  which 
places  certain  constraints  on  material  constants  in  order 
that  the  material  be  admitted  within  the  stress  function. 
Nonetheless,  a large  number  of  materials  satisfy  the 
constraints,  and  Rizzo  and  Shippy  obtained  excellent 
agreement  between  the  numerical  values  obtained  from  a 
piecewise  constant  BIE  (employing  numerical  quadrature 
of  the  necessary  integrals)  and  those  obtained  from 
analytical  solutions  in  direct  comparisons. 

The  BIE  method  has  also  been  applied  to  two 
dimensional  plane  stress  problems  for  fully  anisotropic 
elastic  materials  by  Cruse  [14].  His  development  of  the 
boundary  integral  equations  follows  from  the  notation  and 
theoretical  development  of  the  field  equations  as  given 
by  Lekhnitskii  [15].  Cruse  attained  excellent  agreement 
between  available  exact  analytical  solutions  and  numerical 
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solutions  from  a piecewise  linear  BIE  (utilizing  exact 
integration  of  the  resulting  necessary  integrals) . Note 
that  the  development  of  the  BIE  based  on  the  theoretical 
development  of  Lekhnitskii  requires  that  the  principal 
axes  of  the  material  be  located  in  the  x1~x2  plane  of  the 
problem,  whereas  the  algorithm  of  Clements  and  Rizzo 
allows  the  material  to  assume  a completely  arbitrary 
orientation.  This  arbitrary  orientation  of  the  material 
induces  coupling  of  the  plane  and  antiplane  problems,  so 
that  the  problem  must  be  solved  with  N=3.  Consequently, 
the  size  of  the  coefficient  matrix  (cf.  section  2.3.4) 
necessarily  expands  from  [2n  x 2n]  to  [3n  x 3n] , where  n 
is  the  total  number  of  nodes  in  the  problem's  modelling 
for  the  numerical  solution  by  the  BIE. 

The  plane  stress  problems  of  practical  interest 
usually  involve  the  principal  axes  being  located  in  the 
xl”x2  Plane'  niaking  the  arbitrary  material  orientation 
feature  of  the  algorithm  an  available  convenience  when 
necessary.  A more  practical  application  of  the  case  of 
N=3  is  the  bending  of  plates  weakened  by  openings  under 
antiplane  loading  distributed  along  the  edge.  Lekhnitskii 
(61  gives  extensive  treatment  of  this  type  of  problem, 
which  was  not  considered  in  the  solved  problems  of  this 
thesis. 

One  type  of  plane  stress  problem  of  practical 
interest  is  that  of  plates  weakened  by  openings  subject 
to  various  in-plane  loading  conditions.  For  this  type 
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of  problem  the  variable  of  interest  is  the  hoop  stress 
distribution  around  the  opening.  Inherently  the  boundary 
solution  stress  variables  in  the  BIE  method  are  the 
tractions  from  which  the  hoop  stresses  cannot  be  determined. 
Moreover,  we  cannot  obtain  an  expression  for  3u^/3Xj 
(i=l,2,...N;  j=l, 2)  directly  because  the  numerical 
processes  of  the  PSFBIE  deal  with  the  mapped  coordinate  £. 
However,  we  may  establish  a system  of  simultaneous  linear 
algebraic  equations  to  solve  for  the  stress  system  at 
the  point  in  question.  The  stress  system  may  then  be 
rotated  through  the  necessary  angles  via  Mohr's  circle 
to  obtain  the  oblique  stresses  of  interest.  Appendix  III 
outlines  the  procedure  in  tie  general  case  for  N=2  or  3. 

Finally,  it  should  be  observed  that  the  work  in 
this  thesis  has  been  performed  with  the  idea  that  there 
is  no  absolute  advantage  of  any  numerical  solution  procedure. 
Therefore,  the  need  to  improve  the  approximations  within 
the  numerical  solution  procedures  is  frequently  necessary 
to  obtain  the  desired  solution.  The  foundations  of  these 
improved  approximations  for  the  BIE  method  are  laid  within 
this  thesis  from  which  they  may  be  implemented  readily. 

Also,  it  is  possible  to  generalize  the  methods  presented 
in  this  thesis  to  boundary  value  problems  governed  by 
other  differential  equations,  e.g.,  differential  equations 
which  would  admit  a non-zero  right  hand  side  to  equation 
(2.1)  for  the  case  of  N*1  and  equations  which  govern 
inhomogenous  media  (inclusion  problems)  where  conditions  of 
continuity  at  the  interface  must  be  considered. 

L ' 


' 
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APPENDIX  I 

LIMIT  PROCESSES  FOR  ESTABLISHING  INTERIOR 
IDENTITY  AND  BOUNDARY  INTEGRAL  EQUATION 

I . 1 INTERIOR  IDENTITY 

We  wish  to  establish  an  interior  identity  from 
(2.28)  by  deleting  the  neighborhood  of  the  singular  point 
xq  and  taking  appropriate  limits.  To  facilitate  this 
process,  let  us  assign  explicit  point  dependence  to  the 
variables  involved  in  (2.28). 

► 

[pi(X)*ij  (X,xo)  - rij  (X,xQ)<J)i(X)  ]dS(X)  = 0 (1.1) 


In  classical  fashion  we  delete  a square  region  of 
sides  2e  surrounding  xQ  (see  Figure  1)  whose  contour  we 
will  denote  C' . We  can  rewrite  (1.1)  as 


[pi(X)*ij (X,xo)  - rij(XrxQ)^i(X)]dS(X) 


tPi(5)*ij(5'*o)  ' rij(5'5o)*i(5)lds(?)  " 0 d-2) 


The  first  integral  is  identically  zero  by  our  stipulation 
that  it  does  not  contain  xQ.  Hence  we  need  only  consider 
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the  integral  over  C' . For  convenience  we  let  xQ  be  the 
origin  and  consider  the  parts  of  the  integral  separately. 
Substitution  of  equation  (2,18)  into  (1.2)  with  appropriate 
association  as  outlined  in  section  2.2.1  leads  the  first 
part  to  take  the  form 

pi<?>(5  AiaD0jlo9<2a-°a> 

- £ Aic,5<ijlo9(ic,-5a))dxldx2 


where  the  j subscript  on  Da  has  been  added  to  demonstrate 
the  associations  made  in  section  2.2.1.  Performing  the 
integrations  leads  to  the  expression 


H 


(x|  - 


X.  -c 
1 a 


)log(za-ca)  - 


x -c 
_1 a 

Pa 


2(x  -°a> 


-2  A 


xa  aj 


- I (za-ca)  - - \ 

5 L(X1_C“)  \ ? - x2  j 


(X2  - 


— . 2 


X.  -C 

1 a 


) log (z  -c  ) 
a a 


xrc°V  pcx2_ 

Pc,  /U<x1-ca) 


Substitution  of  the  limits  gives  zero  identically. 
A similar  analysis  on  the  second  part  of  the  integral 


leads  to 


e e 

' _ 
FT  [*i 


iqkl I ^ z -c 
^ \ a a a 


2^i)n 

a z-c  / q 
a a 


+ a.  , _ 
iqk2 


pa  - 2^21  p )„1  ax  dX; 
“ “ “ “ z<j-cc  /•  ^1 


Performing  the  indicated  integrations  leads  to 


the  expression 


i#T  +i<*>  [ aiqkl[?Aic.Daj(  1 °paC‘>iog(za-cc.)  - x2 

tt»  — rr  I ^Za-Ca^  — — \ 

fAkc.Doj( =—■ l-°<3<za-ca>  - x2  J J nq 


+ aiqk2 


[?AicDc 


(Z  -C  ) 

a a 


“D  P r 


l°g(za-ca)  - x2 


__  _ _ /( za-ca)  _ _ r 

|AkaDajPa^ = loglz^cj  - x2jj  nq 


-e  -e 


Substituting  the  designated  limits  with  appropriate  nq 

(i.e. , n = +1  on,  x * +e)  and  letting  e-*-0  gives 

’*1 

the  result  which  contains  the  constant  of  the  integration 


(F)  as 
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*jU0)F 

whereby  (1.1)  assumes  the  desired  form  of  the  interior 
identity 


[P.$ . . - r.  .<p.  ]dS  = -F<f>.  (x  ) 
1 1]  1]  1 j ~o 


I . 2 BOUNDARY  FORMULA 

To  establish  the  boundary  formula,  we  proceed  in 
a similar  fashion  as  in  the  interior  identity  development, 
except  that  rather  delete  the  neighborhood  of  xQ  we 
augment  the  contour.  We  start  with  the  interior  identity 
(2.29)  showing  explicit  point  dependence 

* 

[P. (X)*. . (X,x  ) - r..(X,x  )(fr. (X)]dS(X) 

1 - ij  ~ -O  1]  ~ ~o  1 ~ 

c 

= — F4> j (xQ)  (1.3) 


Augment  boundary  by  part  of  a rectangle  whose 
contour  we  denote  C"  as  shown  in  Figure  2.  Equation  (1.3) 
is  now  realized  by 


“F* j (V 


[P. (X)*.  (X,x  ) 
I ! ~ ij  ~ ~o 

C'+C" 


- rij(x,x0)<>i(x)]ds(x) 


d.4) 


86 


Again  as  in  the  case  of  the  interior  identity,  the 
integral  over  c'  is  identically  zero,  leaving  only  the 
integral  over  C"  to  be  determined.  Closer  examination 
of  the  contour  C"  for  the  case  when  Xq  is  located  at  a point 
on  the  boundary  which  does  not  possess  a unique  tangent 
provides  geometrical  relations  that  will  be  used  in  the 
establishment  of  the  boundary  formula.  Figure  27  depicts 
this  case  for  e small  enough  that  the  neighboring  boundary 
C'  on  either  side  of  Xq  can  be  replaced  with  straight  line 
segments . In  accordance  with  contour  integral  direction 
conventions,  we  denote  the  outward  unit  normal  of  the 
adjacent  straight  line  segment  to  XQ  in  the  positive 
contour  direction  sense  as  n+  and  the  outward  unit  normal 
of  the  adjacent  straight  line  segment  in  the  negative 
contour  direction  sense  as  n~.  We  now  establish  the 
auxiliarly  cartesian  coordinate  system  x'-x'  aligned  as 
shown  in  Figure  27.  Specifically,  x'  is  in  the  -nQ 
direction  where 


n+  + n 


and  x£  is  coincidental  with  the  direction  that  the  contour 
integration  will  follow  along  x'=0. 
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2c 


?o 


Employing  this  auxiliary  cartesian  coordinate  system,  the 
contour  integral  over  C"  becomes 

0 e 


-F*.(xo)  - 


c 

c 


[p.*. . - r.  .<(>.] dxr  + 

i l]  l]  l 2 


[P . $ . . - r . . <J» . ]dx." 
1 1]  1]  1 1 


-e 


(1.5) 


where  z'  = e(l  + tani-(ir-B  (XQ) ) ) . In  the  first  integral 
x£=-e,  n£*-l,  n£=0 , while  in  the  second  integral  x£=0, 
n£=0,  n£=-l,  and  in  the  third  integral  x£=e,  n£=l  and 
n£=0.  We  consider  parts  of  the  integral  separately  so 
that  the  first  parts  of  the  integrals  in  (1.5)  over  C" 
become  by  the  above 
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1 

2iri 


-e 


Pi(?){  l AiaDajlog(_e  + pax2  “ ca) 

# 

- | AiaDajlog(-e  + paX2  - ca) }dx2 
Pi(X){  £ AiaDaj^°9 (xl  “ ca^ 

- £ AiaDajlog(x£  - ca) }dx£ 
pi(X){  E AiaDajlog(e  + paX2  - ca) 

' l A?.aDajlog{£  + pax2  “ ca)  >dx2 


Performing  the  indicated  integrations  gives 
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2iri 
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iaDaj 
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rl  -e+p  x'-c 
— — 1/  2 a 


£ AiaDa] 


log(-e+paX2-ca)-X2 

Pa  / 1 J 1 e 


Pi^?^a  AiaDaj  f ^xl"ca^  ^xl"ca^  “ x£^ 
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pi (X) J E 
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Aiapaj 


2 ct  \ ^ . -• 

log(e+pax2-ca)  - x2 

a / 


~ | AiaDaj|_[ 


17  e+Pax2"Ca 


log (e+paX2~ca)  - xf 
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Putting  in  limits  of  integration  and  letting  e+0 
(whereby  e'-*-0  also)  yields  zero  identically. 

A similar  analysis  on  the  second  parts  of  the 
integrals  in  (1.5)  gives  rise  to  the  expression 


r9 


Jr 


2iri 


LL 


r 


a,  D • 

ka  a' 


A,  D 
ka  a" 


D / Vv  / Kg  ai  _ Kg  g^ \ - 

pi(C)|_aiqkl^  r-e+pax|-ca)  ~ \ (-e+pax£-ca) jnq 


+ a 


iqk2 


Y AkaDajPg AkaDajPg  ^n-ldx^ 

a (~e+Pax2"ca)  a (-£+Pax2"5'a >/  2 


>!<?>[■ 


Y AkaDaj  _ YAkaDctj  \ , 
lrqkl^  a (x£-ca)  a (x£-ca)/  q 


-e 


- 1 


A,  D , 


+ aiqk2 


Aka° 


A,  D . 
ka  aj 


>i(?)  [aiqkl(f  (e+PaX2-ca)  £ (£+Pax2~c^ / M 

(y  AkaDaipa  y Ak°tDctjPct  V. 

aiqk2[2  (E+paX2-Cci)  " i (£+Pa X2-Ca>/  qJ  2 


Performing  the  indicated  integrals  gives 
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f r . A . 

-i-jpitX)  aiqkll^-^^-logt-e+Paxf-Ca) 


- 2 Ct3log(-e+potX2-cct)|  nq 


+ aiqk2  |2AkaDajlog(-e+paX2-ca) 


" 2AkaDajlog(-e+Pcix 


+ Pi<?)[aiqkl(|AkaDaj1°9(xI-ca) 

- 2AkaDajlog(x£-ca)J  n'q 


+ aiqk2l  2AkaDajpalog(x£-ca) 


“ 2AkaDajPalo9(xi-ca))  nq 
a / J 


i(X)  aiqkl  2 -^-^-log ( e+Pax2 -ca ) 
L a Pa 


- Z— =r-*log(e+pax2-ca)lnq 


+ aiqk2(2AkaDajlog(e+pax£-ca) 

^•AkaDaj^°9 ^e+Pax2”ca^|  nqj| 


Substituting  the  designated  limits  with  the 
appropriate  n*  and  letting  e+0  (while  e""e(l+tany(ir-B(X0) ) ) 
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leads  to  the  result 

(BCXq)  - l)F$j(Xo) 

Thus  by  this  result  in  (1.5)  we  arrive  at  the  boundary 
integral  equation 

* 

-Xij(Xo>F*i<x0)  = tPi(x)*i;j(x,x0)  -r^tx.y^widsU) 

c 

where  (XQ)  represents  the  addition  of  the  result  of 
the  contour  integration  over  C"  and  the  left  hand  side 
of  equation  (1.5). 
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APPENDIX  II 


CALCULATION  OF  ELEMENTS  OF  COEFFICIENT  MATRIX 
IN  POLYNOMIAL  SHAPE  FUNCTION  APPROXIMATION 

II. 1 CALCULATION  OF  a“?  AND  FA?j 

We  defined  a^j  and  FA^j  in  section  2.3.4  in 
equations  (2.43)  and  (2.44)  as 


OLO 

a.  . = 
13 


m (5)r.,(x(5) , x )j(5)<U 

i j — -o 


FA? 


rij(X(5),  XQ)J(5)d5 


with  in  the  form  as  given  by  (3.23),  viz. 


rij(X(S),  XQ)  * 2¥I  _^[X!(5)  + PaX2(?)  " (a+pab)] 


2 LikaDaj 

atX^C)  + PaX2(C)  - (a+pab) 

Now  IV  j in  this  form  has  a O(^)  singularity  at  the  point 
ca  ~ 2 a ^ which  suggests  that  we  handle  the  necessary 
integrations  on  the  basis  of  whether  the  singularity 
exists  on  the  segment  CQ  over  which  we  are  integrating. 
Accordingly,  we  can  classify  these  integrals  as  either 


of  the  two  cases: 
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1.  X„  is  not  an  element  of  C 

a 

2.  XQ  is  an  element  of  Ca 

(i)  For  the  first  case,  there  is  no  singularity 
on  Ca  and  the  integrals  can  be  evaluated  approximately  by 
the  v-point  Gaussian  quadrature  rule  given  by  (2.38). 

Thus, 

v 

aij 4 ;0)J(5k),UE 

in  which  A£=2. 

(ii)  For  a polynomial  shape  function  approximation, 
the  occurence  of  XQ  as  an  element  of  Ca  is  recognized 

as  the  case  of  XQ  being  one  of  the  local  nodes  (a)  of  Ca. 

Presupposing  that  $i(XQ)  is  known,  the  behavior 
of  the  left  hand  side  of  equation  (2.39)  may  be  observed 
according  to  the  location  of  XQ.  Specifically  we  want 
to  show  that  the  quantity 

♦i<X>  “ ♦i<?o> 

is  0(r)  (0(C)  in  the  mapped  space)  so  that  the  integrand 

of  the  left  hand  side  of  (2.39)  is  non-singular  since  we 

know  T. .(X,  X ) is  0 (— ) (0(i)  in  the  mapped  space  as  was 
i j - ~o  r £ 

observed  in  the  previous  case) . The  polynomial  shape 
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functions  are  of  the  general  form 


Ma(£)  = aa  + baC  + ca£2  + ... 


Denote  the  value  of  5 at  the  local  node  (a)  at  which  XQ 
is  located  as  £(a).  By  equation  (2.33) 


♦ , (X)  - ♦ , (X  ) 
1 - 1 ~o 


(a°  + baC  + ca£2  + . 


• • ) 0>“° 


(aa  + baC(o)  + ca52(a)  + . . . ) <fr“a 


which  is  0(C).  We  may  rewrite  this  quantity  as 


<K.(X)  - ♦.(Xo)  = ♦®°[M°(*)  - Ma(5(a))] 

i ~ i i 


and  modify  equation  (2.41)  as 


y£l4>iCT  tM<X  (^)  - Ma(C(a))]  (a^  - 


FA?j) 


Z P?°b?? 

CT»1  1 


Define 


a??  - [Ma(C)  - Ma(?(a))](a^  - FA?,) 


as  the  non -singular  matrix  element  to  be  evaluated  by 


the  v-point  Gaussian  quadrature  rule. 
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II. 2 CALCULATION  OF 

a a 

The  matrix  element  b^j  was  defined  in  section  2.3.4 
by  equation  (2.45)  as 

Ma(5)«ij(x(C) ,XQ)J(?)d? 
ca 


with 


*ij(X(C),X0)  - 2iT{|AkoDoj109(Xi<5)+PaX2(5>-(a+Pob)) 

‘ A j 1°9 ( < C ) +PaX2 ( C ) - ( »+Pab> ) } 

As  in  the  previous  section,  there  is  a singularity 
present  in  the  kernel,  which  in  this  case  of  *ij(X(£),X0) 
is  log  (£) . Again  we  separate  cases  of  the  above 
integration  on  the  basis  of  whether  XQ  is  an  element  of  Ca. 

(i)  The  singularity  does  not  exist  on  the  segment 
C0  and  b^j  can  be  evaluated  by  the  v -point  Gaussian 
quadrature  rule  as 

b?j  4 ik£iWk{Ma(5k)*ij(x(ck)  ,x0)j(ek>me 


where  A£«2. 

(ii)  The  order  of  the  singularity  is  log(€)  and 
the  integral  cannot  be  evaluated  in  the  normal  sense. 

To  deal  with  the  situation  we  subdivide  the  segment  Cg 
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into  y subintervals  (where  y is  the  degree  of  the  polynomial 

shape  function  in  the  approximation)  and  integrate  by  an 

appropriate  numerical  procedure  over  the  subintervals 

according  to  whether  XQ  is  a member  of  the  subinterval. 

Denote  the  subintervals  t°  which  are  defined  as 

u> 

t°  * [£(u>)  ,£(<*>+!)  ] o)»l,2,...y 

0) 

in  which  £(<d)  is  interpreted  as  the  value  of  £ at  local 
node  number  u and  £(w+l)  the  value  of  £ at  local  node 
number  w+l.  We  now  have 

Ma (£)♦.. (X(£) »X  ) J(£) d£ 

1 j ■»  ~o 

a 

(i) 

Transform  the  Gaussian  abcissae  £ to  be  in  the  interval 

t°  by 
u > 

_ = te(M)+e(m-i)i  - [£(a>+D -£(u)  i 
4 2 

aa 

We  now  may  approximately  evaluate  b^j  as 

where 


9: 


[S(0))+5(<i>-H)  } - [£<m+1) -£((■>)  ] ^ 


A?  = C(a)+1)  - ?(u>) 


> AC 
'k  = T wk 


If  X0€tJ  we  then  replace  w£  by  the  log  weighting  function  [15] 


w'  = log/(Ck  - C(XQ))2 

where  E(X,J  is  the  value  of  £ at  which  X.  is  located. 
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APPENDIX  III 


CALCULATION  OF  BOUNDARY  STRESSES 


The  mapping  of  the  contour  C to  5 space  does  not 
allow  the  approximation  of  3u^/3Xj  directly  because 
we  cannot  obtain  an  expression  for  3£/3Xj  (cf.  section 
2.3.2).  Consequently,  we  must  establish  a secondary 
system  of  linear  equations  to  evaluate  boundary  values 
of  that  cannot  be  determined  from  the  traction  solution. 
Thereafter,  Mohr's  Circle  can  be  used  to  evaluate  oblique 
stresses  such  as  the  hoop  stress.  The  order  of  the  linear 
system  established  is  dependent  on  the  number  of  components 
(N)  in  the  problem  for  which  appropriate  values  of  indices 
will  be  assigned  subsequently. 

For  the  problem  of  anisotropic  elasticity  the 
fundamental  variables  ( <p are  displacements  (u^)  while 
the  flux-type  variables  (P^)  are  tractions  (t^) . An 
analysis  of  stress  gives  the  relation  [17] 


ti  = a.jnj 


i * 1,2, . . .N 
j * 1,2 


which  represents  N equations  in  2N  unknowns  (a^j) • 
However,  the  anlaysis  of  stress  also  provides  that 


aji 


I 
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so  that  there  are  actually  only  2N-1  unknowns. 

The  assumption  of  small  strains  (Eulerian  strain 
tensor)  allows  formulation  of  N additional  equations. 
The  Eulerian  strain  tensor  €•.,  takes  the  form 


The  PSFBIE  provides  that  by  equation  (2.33)  the  pertinent 
boundary  variables  take  the  form 

u±  = Ma(£)u®  i = 1,2, ...N 

Xj  = Ma(5)x“  j = 1,2 

By  the  chain  rule  for  differentiation  the  derivative  of 
the  displacements  with  respect  to  the  mapped  coordinate 
(£)  can  be  expressed  as 

3ui  3ui  3X. 

aif  * Sx~  u 

Now  the  quantities  3uj/3C  and  3Xj/35  can  be  approximated 
by  use  of  the  shape  functions  (M01)  as 


5ui  , 3Ma(Q  u<* 
H H 1 
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Thus  we  have  established  N equations  in  the  2N  unknowns 
3u^/3Xj,  Note  that  analysis  of  strain  provides  that 
but  not  3u^/3Xj  = 3Uj/3x^,  so  that  at  this 
point  we  have  2N  equations  in  4N-1  unknowns. 

Examination  of  the  constitutive  relation  (Hooke's 
Law)  provides  the  additional  2N-1  equations  in  the  proper 
4N-1  unknowns  necessary  to  complete  the  linear  system. 
Hooke's  Law  is 


°ij  * cijki€kZ 


i, k  = 1,2, ...N 

j, *  - If  2 


where  6^  is  the  Euler ian  strain  tensor.  Using  the 
interpretation  of  £k£  given  above  yields 

Suk  „ 3ui 
°ij  = cijkJl  3x^  + 3xk 

which  may  be  recognized  as  the  required  additional  2N-1 
equations  in  the  proper  4N-1  unknowns. 

Solution  of  the  indicated  linear  system  yields 
values  of  the  stresses  at  the  boundary  point  in  question 
which  may  be  rotated  through  the  necessary  angles  to 
determine  the  stresses  of  interest. 


r 
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